### compile the source code ###
### R CMD SHLIB ??.cpp -DR64 ##
###############################

library(SparseM);

### 3-pt quadrature rule #####
#quadPts<- c(-sqrt(.6), 0, sqrt(.6));
#quadWgts<- c(5/9, 8/9, 5/9) 
############3

####### 5-pt quadrature rule ############3
quadPts<- c(-0.906179845938664,
            -0.538469310105683,
             0.000000000000000,
             0.538469310105683,
             0.906179845938664);
quadWgts<- c(0.236926885056189,
        0.478628670499366,
        0.568888888888889,
        0.478628670499366,
        0.236926885056189);

ems_constructor<- function(centroid, offsets=NULL, bandwidth=0, cellwidth, quadPts, quadWgts) {
	maxNonzero<- (1 + 4*((1+floor(bandwidth/cellwidth))*(2+floor(bandwidth/cellwidth)) - floor(bandwidth*(2-sqrt(2))/(2*cellwidth))^2));
	dimMat<- length(centroid)*maxNonzero;
	if (is.null(offsets)) offsets<- rep(1, length(centroid));
  
	dyn.load("smoother.dll")
	tmpD<- .C("smoother", 
		  cx=as.double(Re(centroid)), 
		  cy=as.double(Im(centroid)),
		  ncell=as.integer(length(centroid)), 
		  quadPts=as.double(quadPts),
		  quadWgt=as.double(quadWgts), 
		  nquadPts=as.integer(length(quadPts)),
		  bandwidth=as.double(bandwidth), 
		  cellwidth=as.double(cellwidth),
		  offsets=as.double(offsets), 
		  resultValue=as.double(rep(0, dimMat)),
		  colNonzero=as.integer(rep(0, length(centroid))),
		  rowIndex=as.integer(rep(-1, dimMat)), 
		  nNonzero=as.integer(0)
		);

	K.mat.csr   <- as.matrix.csr(0);
	K.mat.csr@ra<- tmpD$resultValue[1:tmpD$nNonzero];
	K.mat.csr@ia<- cumsum(c(as.integer(1), tmpD$colNonzero));
	K.mat.csr@ja<- as.integer(1+tmpD$rowIndex[1:tmpD$nNonzero]);
	K.mat.csr@dimension<- rep(length(centroid), 2);
	return(K.mat.csr);

}

intKernel<- function(centroid, bandwidth=0, cellwidth, quadPts, quadWgts) {
	ncell<- as.integer(length(centroid));
	nqdpt<- as.integer(length(quadPts)); 
	nevlt<- ncell*(nqdpt^2 );
	nentry<- (ncell*nqdpt)^2; #ncell^2*(nqdpt^2)

	dyn.load("intKernel.dll")
	tmpD<- 	.C("intKernel", 
		      cx=as.double(Re(centroid)), 
		      cy=as.double(Im(centroid)),
		      ncell=ncell, 
		      quadPts=as.double(quadPts),
		      quadWgt=as.double(quadWgts), 
		      nquadPts=nqdpt,
		      bandwidth=as.double(bandwidth), 
		      cellwidth=as.double(cellwidth),
		      intResult=as.double(rep(0, nentry)), 
		      xevalpt=as.double(rep(0, nevlt)),
		      yevalpt=as.double(rep(0, nevlt))
		);
	
	kernel.mat	<- matrix(tmpD$intResult, ncol=25, byrow=T);

	quadWgts.mat	<- t(diag(1/4, 25) %x% as.vector(outer(quadWgts, quadWgts, "*")));
	
	return(list(kernel.mat=kernel.mat, quadWgts.mat=quadWgts.mat));
}

intEstimate.exact<- function(location, centroid, offsets, bandwidth=0, cellwidth, quadPts, quadWgts) {
  dyn.load("smoother_exact.dll")
  tmpD	<- .C("smoother_exact", cx=as.double(Re(centroid)), cy=as.double(Im(centroid)), 
	ncell=as.integer(length(centroid)), quadPts=as.double(quadPts),
	quadWgt=as.double(quadWgts), nquadPts=as.integer(length(quadPts)),
	bandwidth=as.double(bandwidth), cellwidth=as.double(cellwidth),
	offsets=as.double(offsets), Lambda=as.double(rep(0, length(centroid))),
	xobs=as.double(Re(location)), yobs=as.double(Im(location)), nobs=as.integer(length(location)));

#smoother_exact(double *cx, double *cy, intType *ncell, double *quadPts, double *quadWgt,
#	intType *nquadPts, double *bandwidth, double *cellwidth, double *offsets, double *Lambda,
#	double *xobs, double *yobs, intType *nobs)

}

